#include "HY1C_out.h"
/* ============================================================================ */
/* module l1_pci_hdf.c - functions to read Pathfinder Class data for MSL12      */
/*                                                                              */
/* Modified By: S. Walsh, RSMAS, January 2008 to read Pathfinder                */
/*                                                                              */
/* ============================================================================ */

#include "l1_pci_hdf.h"
#include "l12_proto.h"
#include "hdf_utils.h"
#include "rawcal.h"
#include "runcal.h"

#define NAVBANDS  5
#define NVIRBANDS 2
#define NTIRBANDS 3

#define TEMP_REF 273.16

#define NOA7  7
#define NOA9  9
#define NO10 10
#define NO11 11
#define NO12 12
#define NO14 14
#define NO15 15
#define NO16 16
#define NO17 17
#define NO18 18

//#define SDNAME "Channel_"
#define SDNAME ""


static int32 sd_id;
static int32 sd_id_g;
static int32 spix   = 0;
static int32 sscan  = 0;  /* start scan in original granule */

static float lifetime;


/* ----------------------------------------------------------------------------------- */
/* openl1_pci_hdf() - opens a PCI L1B file for reading.                            */
/*                                                                                     */
/* B. Franz, SAIC, February 2003.                                                      */
/* ----------------------------------------------------------------------------------- */
int openl1_pci_hdf(filehandle *file)
{
    int32  npix;
    int32  nscan;

    int32  sds_id;
    int32  rank;
    int32  dims[3];
    int32  type;
    int32  numattr;

    int32  itemp;
    char   cdata[32];

    /* Open the HDF input file */
    sd_id = SDstart(file->name, DFACC_RDONLY);
    if(sd_id == FAIL){
        fHY1C_out(stderr,"-E- %s line %d: SDstart(%s, %d) failed.\n",
        __FILE__,__LINE__,file->name,DFACC_RDONLY);
        return(HDF_FUNCTION_ERROR);
    }

    /* Get pixel and scan dimensions */
    sds_id = SDselect(sd_id,SDnametoindex(sd_id,(SDNAME"1")));    
    if (SDgetinfo(sds_id,NULL,&rank,dims,&type,&numattr) == -1) {
        fHY1C_out(stderr,"-E- %s line %d: error getting dimension info.\n",
        __FILE__,__LINE__);
        return(HDF_FUNCTION_ERROR);
    }
    npix  = dims[1];
    nscan = dims[0];

    /* Open the HDF geolocation input file.*/

    sd_id_g = SDstart(file->geofile, DFACC_RDONLY);
    if(sd_id_g == FAIL){
        HY1C_out("Error opening geolocation file.\n");
        fHY1C_out(stderr,"-E- %s line %d: SDstart(%s, %d) failed.\n",
            __FILE__,__LINE__,file->geofile,DFACC_RDONLY);
        return(HDF_FUNCTION_ERROR);
    }

    file->npix   = npix;
    file->ndets  = 1;
    file->nscan  = nscan;
    file->sd_id  = sd_id;

    return(LIFE_IS_GOOD);
}


/* ----------------------------------------------------------------------------------- */
/* readl1_pci_hdf() - reads 1 line (scan) from a PCI ingested CLASS file, loads l1rec.        */
/*                                                                                     */
/* ----------------------------------------------------------------------------------- */
int readl1_pci_hdf(filehandle *file, int32 scan, l1str *l1rec)
{
    static int   firstCall = 1;
    static float pi = 3.141592654;

    static  int16 *ch1 = NULL;
    static  int16 *ch2 = NULL;
    static  int16 *ch3 = NULL;
    static  int16 *ch4 = NULL;
    static  int16 *ch5 = NULL;

    static float m   [NAVBANDS];
    static float b   [NAVBANDS];

    static float64 lastTAIsec = -1;
    static float64 TAIsec, usec, dsec;
    static int16  year, jday, RelDay, LY, leap, month, day, hh, mm, ss;
    static int32  ms = 0.0;
    //    static float starttime;
    static char starttime[13];
    static double fsol;
    static int16   mside;
    static double mnorm[3];
    static int32  lastframe = -1;

    static int startOfMonth[2][12] = {{0,31,59,90,120,151,181,212,243,273,304,334},
				      {0,31,60,91,121,152,182,213,244,274,305,335}};
   
    static float rabias;
    static float tempbias;
    static float othbias;

    int16 *data  = NULL; 
    int32 npix   = (int32)file->npix;
    int32 nbands = (int32)file->nbands;
    int32 nbandsir = (int32)file->nbandsir;
    int32 frame  = scan;
    int32 detnum = 1;
    int32 iw,ip,ipb;
    int32 dims[3];

    static float *Fobar;

    char *p1 = NULL;
    char scaletype[128];
    char startdate[15];
    char xsatid[5];

    RAW_CAL *rawrec;
    RUN_CAL *runrec;

    int32 status;

    if (firstCall) {

        firstCall = 0;

        if ((ch1 = (int16 *) calloc(npix, sizeof(int16))) == NULL) {
            HY1C_out("-E- %s line %d: Error allocating data space.\n",
                   __FILE__,__LINE__);
            return(1);
        }

        if (getHDFattr(sd_id,"Slope",SDNAME"1",(VOIDP)&m[0]) != 0) {
            HY1C_out("-E- %s line %d: Error reading Channel_1 Slope attribute.\n",
                   __FILE__,__LINE__);
            return(1);
        }

        if (getHDFattr(sd_id,"Intercept",SDNAME"1",(VOIDP)&b[0]) != 0) {
            HY1C_out("-E- %s line %d: Error reading Channel_1 Intercept attribute.\n",
                   __FILE__,__LINE__);
            return(1);
        }          
        if (getHDFattr(sd_id,"Scale_type",SDNAME"1",(VOIDP) &scaletype) != 0) {
            HY1C_out("-E- %s line %d: Error reading Channel_1 Scale_type attribute.\n",
                   __FILE__,__LINE__);
            return(1);
        }          
	if (strcmp(scaletype,"y= Slope * x + Intercept;") != 0) {
            HY1C_out("-E- %s line %d: Channel_1 Scale_type must be linear.\n",
                   __FILE__,__LINE__);
            return(1);
        }          

        if ((ch2 = (int16 *) calloc(npix, sizeof(int16))) == NULL) {
            HY1C_out("-E- %s line %d: Error allocating data space.\n",
                   __FILE__,__LINE__);
            return(1);
        }

        if (getHDFattr(sd_id,"Slope",SDNAME"2",(VOIDP)&m[1]) != 0) {
            HY1C_out("-E- %s line %d: Error reading Channel_2 Slope attribute.\n",
                   __FILE__,__LINE__);
            return(1);
        }

        if (getHDFattr(sd_id,"Intercept",SDNAME"2",(VOIDP)&b[1]) != 0) {
            HY1C_out("-E- %s line %d: Error reading Channel_2 Intercept attribute.\n",
                   __FILE__,__LINE__);
            return(1);
        }          
        if (getHDFattr(sd_id,"Scale_type",SDNAME"2",(VOIDP) &scaletype) != 0) {
            HY1C_out("-E- %s line %d: Error reading Channel_2 Scale_type attribute.\n",
                   __FILE__,__LINE__);
            return(1);
        }          
	if (strcmp(scaletype,"y= Slope * x + Intercept;") != 0) {
            HY1C_out("-E- %s line %d: Channel_2 Scale_type must be linear.\n",
                   __FILE__,__LINE__);
            return(1);
        }          

        if ((ch3 = (int16 *) calloc(npix, sizeof(int16))) == NULL) {
            HY1C_out("-E- %s line %d: Error allocating data space.\n",
                   __FILE__,__LINE__);
            return(1);
        }

        if (getHDFattr(sd_id,"Slope",SDNAME"3",(VOIDP)&m[2]) != 0) {
            HY1C_out("-E- %s line %d: Error reading Channel_3 Slope attribute.\n",
                   __FILE__,__LINE__);
            return(1);
        }

        if (getHDFattr(sd_id,"Intercept",SDNAME"3",(VOIDP)&b[2]) != 0) {
            HY1C_out("-E- %s line %d: Error reading Channel_3 Intercept attribute.\n",
                   __FILE__,__LINE__);
            return(1);
        }          
        if (getHDFattr(sd_id,"Scale_type",SDNAME"3",(VOIDP) &scaletype) != 0) {
            HY1C_out("-E- %s line %d: Error reading Channel_3 Scale_type attribute.\n",
                   __FILE__,__LINE__);
            return(1);
        }          
	if (strcmp(scaletype,"y= Slope * x + Intercept;") != 0) {
            HY1C_out("-E- %s line %d: Channel_3 Scale_type must be linear.\n",
                   __FILE__,__LINE__);
            return(1);
        }          

        if ((ch4 = (int16 *) calloc(npix, sizeof(int16))) == NULL) {
            HY1C_out("-E- %s line %d: Error allocating data space.\n",
                   __FILE__,__LINE__);
            return(1);
        }

        if (getHDFattr(sd_id,"Slope",SDNAME"4",(VOIDP)&m[3]) != 0) {
            HY1C_out("-E- %s line %d: Error reading Channel_4 Slope attribute.\n",
                   __FILE__,__LINE__);
            return(1);
        }

        if (getHDFattr(sd_id,"Intercept",SDNAME"4",(VOIDP)&b[3]) != 0) {
            HY1C_out("-E- %s line %d: Error reading Channel_4 Intercept attribute.\n",
                   __FILE__,__LINE__);
            return(1);
        }          
        if (getHDFattr(sd_id,"Scale_type",SDNAME"4",(VOIDP) &scaletype) != 0) {
            HY1C_out("-E- %s line %d: Error reading Channel_4 Scale_type attribute.\n",
                   __FILE__,__LINE__);
            return(1);
        }          
	if (strcmp(scaletype,"y= Slope * x + Intercept;") != 0) {
            HY1C_out("-E- %s line %d: Channel_4 Scale_type must be linear.\n",
                   __FILE__,__LINE__);
            return(1);
        }          

        if ((ch5 = (int16 *) calloc(npix, sizeof(int16))) == NULL) {
            HY1C_out("-E- %s line %d: Error allocating data space.\n",
                   __FILE__,__LINE__);
            return(1);
        }

        if (getHDFattr(sd_id,"Slope",SDNAME"5",(VOIDP)&m[4]) != 0) {
            HY1C_out("-E- %s line %d: Error reading Channel_5 Slope attribute.\n",
                   __FILE__,__LINE__);
            return(1);
        }

        if (getHDFattr(sd_id,"Intercept",SDNAME"5",(VOIDP)&b[4]) != 0) {
            HY1C_out("-E- %s line %d: Error reading Channel_5 Intercept attribute.\n",
                   __FILE__,__LINE__);
            return(1);
        }          
        if (getHDFattr(sd_id,"Scale_type",SDNAME"5",(VOIDP) &scaletype) != 0) {
            HY1C_out("-E- %s line %d: Error reading Channel_5 Scale_type attribute.\n",
                   __FILE__,__LINE__);
            return(1);
        }          
	if (strcmp(scaletype,"y= Slope * x + Intercept;") != 0) {
            HY1C_out("-E- %s line %d: Channel_5 Scale_type must be linear.\n",
                   __FILE__,__LINE__);
            return(1);
        }          

	if (getHDFattr(sd_id,"start date","",(VOIDP) &startdate) != 0) {
		HY1C_out("-E- %s line %d: Error reading start date attribute.\n",
		       __FILE__,__LINE__);
		return(1);
	}
	/* convert yyyy-mm-dd to year, day */
	sscanf(startdate,"%4hd-%2hd-%2hd", &year, &month, &day);
	if ( ((year % 400) == 0) || (((year % 4) == 0) && ((year % 100) != 0)) )
	    leap = 1;
	else
	    leap = 0;
	jday = startOfMonth[leap][month-1] + day;
	lifetime = (float) (year - 81) + (float) jday / 365.0;

	if (getHDFattr(sd_id,"start time","",(VOIDP) &starttime) != 0) {
		HY1C_out("-E- %s line %d: Error reading start time attribute.\n",
		       __FILE__,__LINE__);
		return(1);
	}
	/* starttime is a float with a value in this form: hhmmss.ss */
//	hh = starttime/10000.;
//	mm = ((int)(starttime/100.) % 100);
//	ss = ((int)(starttime-(100.0 * (mm+ (100.*hh)))) % 100);
//	ms = (int32_t)(1000. * (starttime- (int)starttime));

	sscanf(starttime,"%2hd:%2hd:%2hd", &hh, &mm, &ss);

	if (getHDFattr(sd_id,"Satellite","",(VOIDP) xsatid) != 0) {
		HY1C_out("-E- %s line %d: Error reading Satellite attribute.\n",
		       __FILE__,__LINE__);
		return(1);
	}

	if (strcmp(xsatid,"NOA7") == 0) l1rec->xsatid = NOA7;
	else if (strcmp(xsatid,"NOA9") == 0) l1rec->xsatid = NOA9;
	else if (strcmp(xsatid,"NO10") == 0) l1rec->xsatid = NO10;
	else if (strcmp(xsatid,"NO12") == 0) l1rec->xsatid = NO12;
	else if (strcmp(xsatid,"NO14") == 0) l1rec->xsatid = NO14;
	else if (strcmp(xsatid,"NO15") == 0) l1rec->xsatid = NO15;
	else if (strcmp(xsatid,"NO16") == 0) l1rec->xsatid = NO16;
	else if (strcmp(xsatid,"NO17") == 0) l1rec->xsatid = NO17;
	else if (strcmp(xsatid,"NO18") == 0) l1rec->xsatid = NO18;

	if (getHDFattr(sd_id,"Satellite Name","",(VOIDP) l1rec->satnam) != 0) {
		HY1C_out("-E- %s line %d: Error reading Satellite Name attribute.\n",
		       __FILE__,__LINE__);
		return(1);
	}

	if (getHDFattr(sd_id,"Orbit","",(VOIDP) &l1rec->morbit) != 0) {
		HY1C_out("-E- %s line %d: Error reading Orbit attribute.\n",
		       __FILE__,__LINE__);
		return(1);
	}

	if (getHDFattr(sd_id,"PRTEMP","",(VOIDP) &l1rec->prtemp) != 0) {
		HY1C_out("-E- %s line %d: Error reading PRTEMP attribute.\n",
		       __FILE__,__LINE__);
		return(1);
	}

	if (getDims(sd_id, "PCI Raw Running Calibration",dims) != 0) {
	    HY1C_out(
		"-E- %s line %d: Error reading PCI Raw Running Calibration dim.\n",
		__FILE__,__LINE__);
	    return(1);
	}
	if ((rawrec = (RAW_CAL *) malloc(dims[0])) == NULL) {
	    HY1C_out("-E- %s line %d: Error allocating rawcal data space.\n",
		   __FILE__,__LINE__);
	    return(1);
	}
	status = rdSDS(sd_id, "PCI Raw Running Calibration",0,0,0,0,(VOIDP) rawrec);

	if (getDims(sd_id, "PCI Running Calibration",dims) != 0) {
	    HY1C_out(
		"-E- %s line %d: Error reading PCI Running Calibration dim.\n",
		__FILE__,__LINE__);
	    return(1);
	}
	if ((runrec = (RUN_CAL *) malloc(dims[0])) == NULL) {
	    HY1C_out("-E- %s line %d: Error allocating runcal data space.\n",
		   __FILE__,__LINE__);
	    return(1);
	}
	status = rdSDS(sd_id, "PCI Running Calibration",0,0,0,0,(VOIDP) runrec);

	if (file->input->degc == 1) {
	    /* Convert to Deg C for calc */
	    rabias = TEMP_REF;
	} else {
	    /* Do calc in Deg K */
	    rabias = 0.;
	}

	tempbias = tempbias - rabias;	/* adjust for C or K */
	othbias  = othbias  - rabias;

	if (leap == 0) {
	    LY=365;
	} else {
	    LY=366;
	}

	switch (l1rec->xsatid) {
	    case NOA7:
		HY1C_out("I Applying NOAA-7 visible degradation.\n");
		m[0] = m[0]/(0.916-0.049*(year+jday/LY-81.5)+
			0.0050*pow(year+jday/LY-81.5,2));
		b[0] = b[0]/(0.916-0.049*(year+jday/LY-81.5)+
			0.0050*pow(year+jday/LY-81.5,2));
		m[1] = m[1]/(0.882-0.080*(year+jday/LY-81.5)+
			0.0110*pow(year+jday/LY-81.5,2));
		b[1] = b[1]/(0.882-0.080*(year+jday/LY-81.5)+
			0.0110*pow(year+jday/LY-81.5,2));
		break;
	    case NOA9:
		HY1C_out("I Applying NOAA-9 visible degradation.\n");
		m[0] = m[0]/(0.953-0.051*(year+jday/LY-85));
		b[0] = b[0]/(0.953-0.051*(year+jday/LY-85));
		m[1] = m[1]/(0.866-0.026*(year+jday/LY-85));
		b[1] = b[1]/(0.866-0.026*(year+jday/LY-85));
		break;
	    case NO11:
		HY1C_out("I Applying NOAA-11 visible degradation.\n");
		m[0] = m[0]/(0.797-0.010*(year+jday/LY-89));
		b[0] = b[0]/(0.797-0.010*(year+jday/LY-89));
		m[1] = m[1]/(0.683-0.020*(year+jday/LY-89));
		b[1] = b[1]/(0.683-0.020*(year+jday/LY-89));
		break;
	    case NO14:
		/* Dec 30, 1994 is base date */
	        RelDay = ((year-94)*365+(jday-364));
		HY1C_out("I Applying NOAA-14 visible degradation.\n");
		/* Post-Launch Calibration for Noaa-14 */
		/* C R Nagaraja Rao and Jianhua Chen   */
		/* NOAA/NESDIS                         */
		/* (from web page, March 1997)         */
		m[0] = 0.0000232*RelDay + 0.109;
		b[0] = -41.0*m[0];
		m[1] = 0.0000373*RelDay + 0.129;
		b[1] = -41.0*m[1];
		break;
	    default:
		HY1C_out("I No visible degradation available for NOAA-%2.0d.\n",
			l1rec->xsatid);
		break;

	}

    }

    *(l1rec->year) = (int32_t)year;
    *(l1rec->day)  = (int32_t)jday;
    *(l1rec->msec) = (int32_t)(((hh*60+mm)*60+ss)*1000+ms);

    /* Get position and path geometry */
    READ_SDS_ID(sd_id_g,"lon",l1rec->lon, scan,spix,0,0, 1,npix,1,1);
    READ_SDS_ID(sd_id_g,"lat" ,l1rec->lat, scan,spix,0,0, 1,npix,1,1);
    READ_SDS_ID(sd_id_g,"solz" ,l1rec->solz, scan,spix,0,0, 1,npix,1,1);
    READ_SDS_ID(sd_id_g,"satz" ,l1rec->senz, scan,spix,0,0, 1,npix,1,1);
    READ_SDS_ID(sd_id_g,"relazim",l1rec->delphi, scan,spix,0,0, 1,npix,1,1);

    for (ip=0; ip<npix; ip++) {

        l1rec->pixnum[ip] = spix + ip;

        if (l1rec->lon[ip] < -181.0 || l1rec->lon[ip] > 181.0 || 
            l1rec->lat[ip] <  -91.0 || l1rec->lat[ip] >  91.0 )
	  l1rec->navfail[ip] = 1; 
    }


    /* Read L1B data, scale to radiance, and copy relevant bands to L1 record  */
    READ_SDS_ID(sd_id,SDNAME"1",ch1, scan,spix,0,0, 1,npix,1,1);
    READ_SDS_ID(sd_id,SDNAME"2",ch2, scan,spix,0,0, 1,npix,1,1);
    READ_SDS_ID(sd_id,SDNAME"3",ch3, scan,spix,0,0, 1,npix,1,1);
    READ_SDS_ID(sd_id,SDNAME"4",ch4, scan,spix,0,0, 1,npix,1,1);
    READ_SDS_ID(sd_id,SDNAME"5",ch5, scan,spix,0,0, 1,npix,1,1);

    for (iw=0; iw<nbands; iw++) {

        switch (iw) {
          case 0: data = ch1; break;
          case 1: data = ch2; break;
        }

        for (ip=0;ip<npix; ip++) {

	    ipb = ip*NBANDS+iw;
	    l1rec->Lt[ipb] = 0.0;

	    /* check for sentinel values and flag as appropriate */
  	    if (data[ip] >= 32000) {

  	        switch (data[ip]) {

		default:
		  l1rec->hilt[ip] = 1;
		  l1rec->Lt[ipb] = 1000.0;
		  break;
	        }
		
	    } else {
	        l1rec->Lt[ipb] = (data[ip] - b[iw]) * m[iw];
	    }
  	}
    }

    for (iw=0; iw<nbandsir; iw++) {

        switch (iw) {
          case 0: data = ch3; break;
          case 1: data = ch4; break;
          case 2: data = ch5; break;
        }

        for (ip=0;ip<npix; ip++) {

	    ipb = ip*NBANDS+iw;
	    l1rec->Lt[ipb] = 0.0;

	    /* check for sentinel values and flag as appropriate */
  	    if (data[ip] >= 32000) {

  	        switch (data[ip]) {

		default:
		  l1rec->hilt[ip] = 1;
		  l1rec->Lt[ipb] = 1000.0;
		  break;
	        }
		
	    } else {
	        l1rec->Ltir[ipb] = (data[ip] - b[iw]) * m[iw];
	    }
  	}
    }

    l1rec->sensorID = file->sensorID;
    l1rec->npix     = file->npix;
    l1rec->detnum   = (int32_t)detnum;
    l1rec->mside    = (int32_t)mside;

    /* Convert IR bands to brightness temperature */
    //    radiance2bt(l1rec,-1);    

    return(LIFE_IS_GOOD);
}


int closel1_pci_hdf(filehandle *file)
{
   if ((sd_id_g != sd_id) && SDend(sd_id_g)) {
      fHY1C_out(stderr,"-E- %s line %d: SDend(%d) failed for file, %s.\n",
      __FILE__,__LINE__,sd_id,file->geofile);
      return(HDF_FUNCTION_ERROR);
   }
   if (SDend(sd_id)) {
      fHY1C_out(stderr,"-E- %s line %d: SDend(%d) failed for file, %s.\n",
      __FILE__,__LINE__,sd_id,file->name);
      return(HDF_FUNCTION_ERROR);
   }

   return(LIFE_IS_GOOD);
}
 
 

